Kauffman Boolean model in undirected scale free networks 
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We investigate analytically and numerically the critical line in undirected random Boolean net- 
works with arbitrary degree distributions, including scale- free topology of connections P{k) ^ fe~^. 
We show that in infinite scale-free networks the transition between frozen and chaotic phase occurs 
for 3 < 7 < 3.5. The observation is interesting for two reasons. First, since most of critical phenom- 
ena in scale-free networks reveal their non-trivial character for 7 < 3, the position of the critical line 
in Kauffman model seems to be an important exception from the rule. Second, since gene regula- 
tory networks are characterized by scale-free topology with 7 < 3, the observation that in finite-size 
networks the mentioned transition moves towards smaller 7 is an argument for Kauffman model as 
a good starting point to model real systems. We also explain that the unattainability of the critical 
line in numerical simulations of classical random graphs is due to percolation phenomena. 
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Almost 40 years ago Stuart Kauffman proposed ran- 
dom Boolean networks (RBNs) for modelling gene regu- 
latory networks [ij. Since then, beside its original pur- 
pose, the model and its modifications have been applied 
to many different phenomena like cell differentiation [3 , 
immune response 0] , evolution 0] , opinion formation , 
neural networks 6] , and even quantum gravity problems 
0. 

The original RBNs were represented by a set of N el- 
ements, = {cri(t), (T2(i), crjv(t)}, each element Ui 
having two possible states: active (1), or inactive (0). 
The value of at was controlled by k other elements of 
the network, i.e. Uiit + 1) = fi{ai-,{t),CTi^{t), ...,i7i^{t)), 
where k was a fixed parameter. The functions fi were 
selected so that they have returned values 1 and with 
probabilities respectively equal to p and 1 — p. The pa- 
rameters k and p have determined the dynamics of the 
system (Kauffman network) , and it has been shown that 
for a given probability p, there exists the critical number 
of inputs [13| 
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below which all perturbations in the initial state of the 
system die out (frozen phase), and above which a small 
perturbation in the initial state of the system may prop- 
agate across the entire network (chaotic phase). 

In fact, the behavior of Kauffman model in the vicin- 
ity of the critical line kc(p) has become a major con- 
cern of scientists interested in gene regulatory networks. 
The main reason for this was the conjecture that living 
organisms operate in a region between order and com- 
plete randomness or chaos (the so-called edge of chaos) 
where both complexity and rate of evolution are maxi- 
mized [1, . The analogous behavior has been noticed 
in Kauffman networks, which in the interesting region de- 
scribed by eq. ^ show stability, homcostatis, and the 
ability to cope with minor modifications when mutated. 
The networks are stable as well as flexible in this region. 



Recently, when data from real networks have become 
available [111 , [l^ , a quantitative comparison of the edge 
of chaos in these datasets and RBN models has brought 
an encouraging and promising message that even such 
simple model may quite well mimic characteristics of real 
systems. 

Since, however, one has noticed that real genetic net- 
works exhibit a wide range of connectivities, the recent 
modifications of the standard RBN take into considera- 
tion a distribution of nodes' degrees P(k). It has been 
shown that if the random topology of the directed net- 
work is homogeneous (i.e. all elements of the network are 
statistically equivalent), then the network topology can 
be meaningfully characterized by the average in-degree 
(fc) , and the transition between frozen and chaotic phase 
occurs for [1411: 
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On the other hand, if the network topology is char- 
acterized by a wide heterogeneity in the connectivity of 
elements, then it is useless to characterize the network by 
the average in-degree, and instead of (fc) another param- 
eter must be used. In the case of power-law in-degree 
distribution P(k) = [C('^)k'']^'^ , where ^(7) is the zeta 
function, the characteristic exponent 7 is the relevant 
parameter. It has been shown that the critical line 7c (p) 
in RBN model defined on scale-free networks is given by 



C(7c 



1) 
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C(7c) 2p(l-p)- 



(3) 



Since 2 < jc{p) < 2.5, based on the result ([3]) it was 
claimed 16] that the abundance of scale-free networks 
with 2 < 7 < 3 in nature and society can be attributed 
to the presence of both phases, frozen and chaotic, only 
is such networks. 

Recently, several authors 0, 0] have provided a gen- 
eral formula for the edge of chaos in directed networks 
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characterized by the joint degree distribution P{k, q) 
(kq) 1 



{q) 2p(l-p)' 



(4) 



where k and q correspond to in- and out-degrees of the 
same node, respectively. The formula ([4]) shows that the 
position of the critical line depends on the correlations 
between k and q in such networks. It is also easy to 
show that the previous results ([I])-® immediately follow 
from ([3]) if one assumes the lack of correlations P{k, q) — 

P^n{k)Pout{q)■ 

In this paper, we derive general relation describing po- 
sition of the critical line in undirected RBNs with ar- 
bitrary distribution of connections P{k). The specific 
cases, including homogeneous as well as strongly hetero- 
geneous (i.e. scale-free) random network topologies are 
discussed. We also generalize our derivations to the case 
when the scale- free network topology is characterized not 
only by the exponent 7 but also by the minimal node de- 
gree kmin = TTL, which controls the density of connections. 
We show that for 7 — > 00 the parameter m corresponds 
to the original parameter k used in the standard Kauff- 
man model defined on regular random graphs, in which 
the number of connections is the same for all elements. 

In order to find the position of the critical line in RBN 
one has to examine the sensitivity of its dynamics with re- 
gard to the initial conditions. In numerical studies such a 
sensitivity can be analyzed quite simply. One has to start 
with two initial states — {cri(O), (72(0), crAr(O)} and 

X^o ~ {2^1(0), ^2(0), ...,crAr(0)}, which are identical ex- 
cept for a small number of elements, and observe how 
the differences between both configurations and 
change in time. If a system is robust then the studied 
configurations lead to similar long-time behavior, other- 
wise the differences develop in time. A suitable measure 
for the distance between the configurations is the overlap 
x{t) defined as 



x{t) 



1 ^ 



(5) 
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Note, that in the limit N 00, the overlap becomes the 
probability for two arbitrary but corresponding elements, 
ai{t) and a^i{t), to be equal. Moreover, the stationary 
long-time limit of the overlap x — limt^oo x(t) can be 
treated as the order parameter of the system. If a; = 
1 then the system is insensitive to initial perturbations 
(frozen phase), while for x < 1, the initial perturbations 
propagate across the entire network (chaotic phase) . 

In the following, we will partially reproduce the an- 
nealed computation (for the first time carried out by Der- 
rida and Pomeau pJi]), and generalize it to the case of 
undirected random graphs with arbitrary degree distri- 
bution. The case of directed networks has been studied 
by Aldana , and also by Lee and Rieger [l3| ■ 

Thus, having in mind that x(t) corresponds to the 
probability that a given element i possesses the same 



value in both configurations, <7i(t) = ^i{t), two different 
situations have to be considered. If all the ki inputs of 
ai{t) are equal to respective inputs oiai{t), which occurs 
with probability [^(i)]'^', then one has cr,i(t-|-l ) = ai{t+l). 
On the other hand, if at least one of the ki inputs of crj(t) 

differs from its counterpart in which occurs with 

probability 1 — [a;(t)]'^', then ai(t + 1) = ai{t + 1) only 
if /i(CTji(t),...,crjfc. (t)) = fi{ai,{t),...,ai^^{t)) regardless 
of the values of the inputs in each configuration. Prob- 
ability of such an event is + {1 — . Taking all the 
above together one finds that the probability x{t^l) that 
ai{t + 1) = cri(< + 1) is given by 
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{\-[x{tt^)-{p^^{\-pf)}Q{h), (6) 



where Q{ki) represents probability that an arbitrary link 
leads to the node i of degree ki. In uncorrelated networks 
Q{k) corresponds to the degree distribution of the nearest 
neighbors 



(7) 



The equation ([6]) can be understood as a map x{t+l) = 
M{x{t)), where 

00 

M{x) = l-2p{l-p){l-J2x''Q{k)} (8) 

fc=i 

It can be shown that the change of the stability of the 
fixed point of the map x — M{x), which occurs when 



dM{x) 
dx 



(9) 



determines the phase transition between the ordered and 
chaotic regimes (c.f. [l^)- Substituting ([5]) into ^ one 
gets the condition for the phase transition: 
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(10) 



In the following we will analyze the equation (|10p in clas- 
sical random graphs and in scale- free networks where the 
second moment (fc^) becomes important (it diverges for 
7<3). 

Since in classical random graphs (fc^) = (fc)^ -I- (fc), the 
eq. (fTU|) simplifies: 
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Comparing the formula with ^ one can see that 
the critical curve in undirected networks has been shifted 
by 1 in comparison with the directed case. The figure [1] 
presents both equations as well as numerical simulations 
of undirected networks of three different sizes [N = 50, 
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FIG. 1: Phase diagram of Kauffman model defined on clas- 
sical random graphs. Frozen phase resulting from eq. (|lip 
is marked by the dark gray area. Light gray area shows the 
difference between directed and undirected networks. 



N = 500, and N = 5000). While in the limit of large 
{k) the results, especially for large N, agree very well 
with the eq. ([TT|) (see inset), for (k) 1 (i.e. p 0.5) 
they differ significantly. The discrepancy results from the 
fact that (k) = 1 corresponds to the percolation thresh- 
old in these networks. Because the size of the largest 
component near (fc) = 1 is significantly smaller than the 
network size (the network is divided into several not con- 
nected components) , any perturbation cannot propagate 
across the entire system, and the frozen phase is easier 
achieved. It means that it is impossible to verify eq. (fTTj) 
in this range. The closer percolation threshold we are the 
smaller networks (separated pieces of the whole network) 
we analyze. One can also show that if one introduces as- 
sortativity (i.e. positive degree-degree correlations) to 
the network the attainable critical connectivity can be 
significantly shifted towards (fc)c = 1. It happens be- 
cause the percolation transition occurs for lower values 
of (k) in assortative networks [13, [2l|. Unfortunately, 
due to the introduced correlations, analytical treatment 
is much more difficult in such a case. 

Now, let us analyze scale- free networks with the degree 
distribution given by power law 



P(fc) = [C(7,m)F] 



(12) 



where ^(7, m) = Sfe^m generalized Riemann 

zeta function (normalization factor), and the parameter 
m represents the minimal node degree, i.e. it controls 
density of connections in the considered networks. Now 
the eq. PU)) takes a form: 



(13) 



C(7-2,m) ^ 1 
C(7-l,m) 2p{l-p)- 



In figured] comparison of transcendental equations ([T 
(undirected network) for m = 1 and Q (directed net 



FIG. 2: Phase diagram of scale- free networks with m — 1. 
Frozen phase resulting from eq. (I13|l is marked by dark gray 
area. Light gray area shows the difference between directed 
and undirected networks. Points represent results of numeri- 
cal simulations, while the two intermediate lines are solutions 
of eq. p3|l modified for finite networks (dot-dashed line for 
N = 500, and dotted line for N = 5000). 
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FIG. 3: Critical lines for scale-free RBNs with m = 1 (solid 
lines and filled points) , m — 2 (dashed lines and open points) 
and m = 3 (dotted lines and crosses). Thick lines are solutions 
of eq. (I13|l . while the thin lines are solutions of the same 
equation modified for networks of size A'' = 5000. Points 
correspond to the results of numerical simulations. 



work) is presented. Analytical curves taking into account 
finite size version of the distribution (|12p (where zeta 
functions have been replaced by finite sums), as well as 
results of the numerical simulations for N — 500 and 
N — 5000 are also shown in the figure. One can see 
that in undirected case of infinite scale-free networks the 
transition between frozen and chaotic phase occurs for 
3 < 7 < 3.5. It means that in the studied network the 



Y 



30 
251- 
20 
151- 
10 
51- 




_] I I I I I L_ 



_] I L_ 



0.0 0.2 0.4 0.6 0.8 1.0 
P 

FIG. 4: Phase diagram of scale-free networks with different 
values of the parameter m. 



critical line has been shifted in comparison with the di- 
rected case by A7 = 1 towards larger values of the expo- 
nent 7. 

The observation is interesting for two reasons. First, 
since gene regulatory networks are characterized by scale- 
free topology with 7 < 3, the observation that in finite- 
size networks the mentioned transition moves towards 
smaller 7 is an argument for Kauffman model as a good 
starting point to model real systems. Second, most of 
critical phenomena in scale-free networks reveals its non- 
trivial character for 7 < 3 making these networks inter- 
esting for researchers [22] . It happens because the second 
moment of the degree distribution is size dependent for 
7 < 3 (it diverges for N 00). For example, in the 
case of percolation transition, the above causes that it 
is practically impossible to eliminate the giant connected 
component in such networks, i.e. they are ultraresilient 
against random damage or failures [23, [IJ] . It also im- 
plies the lack of epidemic threshold in such networks, i.e. 
the networks are prone to the spreading and the persis- 
tence of infections whatever the epidemic spreading rate 
is. Finally, in Ising model defined on scale-free networks 
with 7 < 3 the critical temperature is size dependent. 
Taking all the above into consideration the position of 
the critical line in Kauffman model shows that scale-free 
networks with 7 > 3 may also exhibit interesting prop- 
erties. 

In previous papers [l^ [lB| it has been stated that the 
only natural parameter which determines the network 
topology is the scale-free exponent 7. In this paper, we 
introduce the parameter m, which does not change the 
scale-free character of the node degree distribution, but 
allows us to control the density of connections. For m — 1 
we retrieve the original problem studied in [isl Il6| . In 
figure [3] and |4] we present the solutions of the eq. (|13p for 
different values of the parameter m. As one can see, for 
m > 2 the frozen phase is preserved only for sufficiently 
small and for sufficiently high values of the parameter p. 
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FIG. 5: Phase diagram of scale-free RBNs with 7 00. The 
diagram coincides with the phase diagram of the standard 
RBN model. Points for a given m show the width of the 
chaotic phase taken from the fig. |3]for 7 = 30. 



For a wide range of intermediate values of p the frozen 
phase is unattainable. 

It is worth noting that in the limit 7 — > 00, the scale- 
free distribution (|12p transforms into the Dirac delta 
function S{k — m) (then (fc^) = (fc)^ and eq. pn]) sim- 
plifies to eq. ([!])). It means that in this limit the scale- 
free RBN model transforms to the standard RBN model, 
where all elements have the same node degree. In fig. |3] 
one can see that for 7 — > cjo and m = 2 the width of the 
chaotic phase shrinks to zero. In figure [5] we show this 
width for different values of the parameter m. In this 
figure one can easily recognize the phase diagram of the 
standard RBN model, in which for p ~ 0.5 the critical 
value of the node degree equals to to = fcc = 2. 

In summary, we have investigated analytically and nu- 
merically the critical line in undirected random Boolean 
networks with arbitrary degree distribution including ho- 
mogeneous and scale-free topology of connections. We 
have shown that in infinite scale- free networks the transi- 
tion between frozen and chaotic phase occurs for 3 < 7 < 
3.5, i.e. position of the critical line is shifted by A7 — 1 
towards larger values of the exponent 7 in comparison 
with the directed case. The observation is interesting 
for two reasons. First, since most of critical phenomena 
in scale-free networks reveals its non-trivial character for 
7 < 3, the position of critical line in Kauffman model 
seems to be an important exception from the rule. Sec- 
ond, since gene regulatory networks are characterized by 
scale- free topology with 7 < 3, the observation that in 
finite-size networks the mentioned transition moves to- 
wards smaller 7 is an argument for Kauffman model as 
a good starting point to model real systems. We also 
explain that the unattainability of the critical line in nu- 
merical simulations of classical random graphs is due to 
percolation phenomena. 
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